Measurement exploratory data analysis

DEVICES = ['icomox', 'icomox-fast', 'beaglebone', 'beaglebone-bvs', 'beaglebone-old', 'beaglebone-fan']

DEVICE = DEVICES[3]
T_WAVEFORM = 5
T_SEC = T_WAVEFORM
NFFT = 1024

from tabulate import tabulate
from IPython.display import Markdown, HTML
import matplotlib.pylab as plt
from typing import List, Tuple

import pandas as pd
import numpy as np
from scipy.signal import find_peaks
from scipy.signal import butter, lfilter

from tsfel.feature_extraction.features import fundamental_frequency

import os
import os.path
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda
from vibrodiagnostics import discovery
from vibrodiagnostics import datasets

pd.set_option("display.notebook_repr_html", False)

file_list = []

if DEVICE == 'icomox':
    Fs = 400
    DIR = '../../air_conditioning/shc3/'
    file_list = ['icomox_k5_low_power_accel.txt']
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.icomox_dataset(file_list)

elif DEVICE == 'icomox-fast':
    Fs = 50000
    DIR = '../../air_conditioning/shc3/'
    file_list = ['icomox_k5_accel.txt']
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.icomox_dataset(file_list)

elif DEVICE == 'beaglebone':
    Fs = 2500
    file_list = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    DIR = '../../air_conditioning/shc3/'
    
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-bvs':
    Fs = 2500
    file_list = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv', 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]
    DIR = '../../air_conditioning/bvs/'
    
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-old':
    Fs = 1160     # 1160 Hz +/-150 Hz (1020 - 1340 Hz)
    file_list = [
        'excool_digitalis/high_pressure_pump_40percent__1.tsv',
        'excool_digitalis/high_pressure_pump_40percent__2.tsv',
        'excool_digitalis/high_pressure_pump_40percent_speed_up__3.tsv',
        'excool_digitalis/high_pressure_pump_80percent__4.tsv',
        'shc2/compressor_shc2_klima4_bad.tsv',
        'shc2/compressor_shc2_klima5_good.tsv',
        'vertiv_digitalis/compressor_top_1.tsv',
        'vertiv_digitalis/compressor_top_2.tsv',
        'vertiv_digitalis/compressor_top_3.tsv',
        'vertiv_digitalis/compressor_base.tsv',
        'vertiv_digitalis/compressor_side.tsv'
    ]
    DIR = '../../air_conditioning/'
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-fan':
    Fs = 2500
    file_list = [
        'fan/speeds/1_still.tsv',
        'fan/speeds/2_still.tsv',
        'fan/speeds/3_still.tsv',
        'fan/speeds/1_up.tsv',
        'fan/speeds/2_up.tsv',
        'fan/speeds/3_up.tsv',
        'fan/speeds/1_down.tsv',
        'fan/speeds/2_down.tsv',
        'fan/speeds/3_down.tsv',
    ]
    DIR = '../../air_conditioning/'
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = beaglebone_dataset(file_list)

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

bvs_1_hore.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 76500 entries, 0.0 to 30.599600000000002
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       76500 non-null  float64
 1   y       76500 non-null  float64
 2   z       76500 non-null  float64
dtypes: float64(3)
memory usage: 2.3 MB

bvs_2_hore.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

bvs_3_motor.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 76500 entries, 0.0 to 30.599600000000002
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       76500 non-null  float64
 1   y       76500 non-null  float64
 2   z       76500 non-null  float64
dtypes: float64(3)
memory usage: 2.3 MB

bvs_4_motor.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))

bvs_1_hore.tsv

	x	y	z
count	76500	76500	76500
mean	4.55518e-15	1.23695e-14	3.69863e-15
std	0.39246	0.6742	0.483534
min	-1.53161	-2.60494	-1.96145
25%	-0.262246	-0.473369	-0.332838
50%	0.00120628	0.00563534	0.026415
75%	0.264658	0.460689	0.337768
max	1.70167	2.54436	1.70293

bvs_2_hore.tsv

	x	y	z
count	51000	51000	51000
mean	6.45201e-16	1.4482e-14	7.69991e-15
std	0.396248	0.674178	0.523105
min	-1.54745	-2.62487	-1.93246
25%	-0.278093	-0.469353	-0.375698
50%	-0.0146411	-0.0142992	0.00750486
75%	0.272761	0.464705	0.366758
max	1.66187	2.73997	2.21092

bvs_3_motor.tsv

	x	y	z
count	76500	76500	76500
mean	-5.27102e-15	6.92227e-15	-4.63692e-15
std	0.174639	0.223881	0.166772
min	-0.674186	-0.822784	-0.541676
25%	-0.123332	-0.176129	-0.110573
50%	-0.00358063	0.0154725	0.00917809
75%	0.11617	0.159174	0.104979
max	0.834676	0.686078	0.655833

bvs_4_motor.tsv

	x	y	z
count	51000	51000	51000
mean	7.84342e-15	8.43428e-15	3.96232e-16
std	0.194656	0.224898	0.165745
min	-0.822007	-0.695444	-0.611263
25%	-0.127452	-0.16854	-0.108309
50%	-0.00770069	0.0230617	0.0114416
75%	0.136	0.166763	0.107242
max	0.806606	0.741568	0.634147

Time domain histogram

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].boxplot(figsize=(8, 3), grid=True)
    plt.show()

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

Time domain waveform

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

Time domain waveform zoom detail

for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

def spectogram(x):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))

Spectrogram in X axis

x_spectra = get_spectrograms(DATASET, 'x')

bvs_1_hore.tsv (X axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_2_hore.tsv (X axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_3_motor.tsv (X axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_4_motor.tsv (X axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in X axis

show_spectrogram_detail(x_spectra, 'x', T_SEC)

bvs_1_hore.tsv (X axis @ 5s)

bvs_2_hore.tsv (X axis @ 5s)

bvs_3_motor.tsv (X axis @ 5s)

bvs_4_motor.tsv (X axis @ 5s)

Peaks in frequency spectrum in X axis

• MMS peak finder algorithm

show_mms_peaks(x_spectra, 'x', T_SEC)

bvs_1_hore.tsv (X axis @ 5s)

bvs_2_hore.tsv (X axis @ 5s)

bvs_3_motor.tsv (X axis @ 5s)

bvs_4_motor.tsv (X axis @ 5s)

Harmonic series detection in X axis

# show_harmonic_series(x_spectra, 'x', T_SEC)

show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

bvs_1_hore.tsv (X axis @ 5s)

• Fundamental frequency: 63.4765625 Hz

	f	y	1x
0	63.4766	0.113014	1
1	85.4492	0.0379399	1.34615
2	129.395	0.0272333	2.03846
3	593.262	0.0227783	9.34615
4	541.992	0.0186733	8.53846

bvs_2_hore.tsv (X axis @ 5s)

• Fundamental frequency: 61.03515625 Hz

	f	y	1x
0	61.0352	0.0873599	1
1	146.484	0.0349362	2.4
2	100.098	0.0218704	1.64
3	183.105	0.0201232	3
4	1098.63	0.0173606	18

bvs_3_motor.tsv (X axis @ 5s)

• Fundamental frequency: 75.68359375 Hz

	f	y	1x
0	75.6836	0.0319051	1
1	24.4141	0.0234566	0.322581
2	292.969	0.0145856	3.87097
3	170.898	0.00893673	2.25806
4	1022.95	0.0074295	13.5161

bvs_4_motor.tsv (X axis @ 5s)

• Fundamental frequency: 97.65625 Hz

	f	y	1x
0	97.6562	0.049372	1
1	107.422	0.0299587	1.1
2	1123.05	0.0203222	11.5
3	134.277	0.0121812	1.375
4	393.066	0.00889755	4.025

Spectrogram in Y axis

y_spectra = get_spectrograms(DATASET, 'y')

bvs_1_hore.tsv (Y axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_2_hore.tsv (Y axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_3_motor.tsv (Y axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_4_motor.tsv (Y axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Y axis

show_spectrogram_detail(y_spectra, 'y', T_SEC)

bvs_1_hore.tsv (Y axis @ 5s)

bvs_2_hore.tsv (Y axis @ 5s)

bvs_3_motor.tsv (Y axis @ 5s)

bvs_4_motor.tsv (Y axis @ 5s)

Peaks in frequency spectrum in Y axis

show_mms_peaks(y_spectra, 'y', T_SEC)

bvs_1_hore.tsv (Y axis @ 5s)

bvs_2_hore.tsv (Y axis @ 5s)

bvs_3_motor.tsv (Y axis @ 5s)

bvs_4_motor.tsv (Y axis @ 5s)

Harmonic series detection in Y axis

# show_harmonic_series(y_spectra, 'y', T_SEC)

show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

bvs_1_hore.tsv (Y axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.252131	1
1	468.75	0.0660276	19.2
2	593.262	0.056911	24.3
3	402.832	0.0511265	16.5
4	603.027	0.0443939	24.7

bvs_2_hore.tsv (Y axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.244484	1
1	583.496	0.0807483	23.9
2	405.273	0.054666	16.6
3	646.973	0.0454914	26.5
4	678.711	0.0403614	27.8

bvs_3_motor.tsv (Y axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.119071	1
1	925.293	0.0126255	37.9
2	231.934	0.00894881	9.5
3	68.3594	0.00792483	2.8
4	131.836	0.0067397	5.4

bvs_4_motor.tsv (Y axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.131803	1
1	56.1523	0.0104745	2.3
2	302.734	0.00907972	12.4
3	283.203	0.00819615	11.6
4	222.168	0.00755379	9.1

Spectrogram in Z axis

z_spectra = get_spectrograms(DATASET, 'z')

bvs_1_hore.tsv (Z axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_2_hore.tsv (Z axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_3_motor.tsv (Z axis, Fundamental = 24.3791 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

bvs_4_motor.tsv (Z axis, Fundamental = 24.3627 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Z axis

show_spectrogram_detail(z_spectra, 'z', T_SEC)

bvs_1_hore.tsv (Z axis @ 5s)

bvs_2_hore.tsv (Z axis @ 5s)

bvs_3_motor.tsv (Z axis @ 5s)

bvs_4_motor.tsv (Z axis @ 5s)

Peaks in frequency spectrum in Z axis

show_mms_peaks(z_spectra, 'z', T_SEC)

bvs_1_hore.tsv (Z axis @ 5s)

bvs_2_hore.tsv (Z axis @ 5s)

bvs_3_motor.tsv (Z axis @ 5s)

bvs_4_motor.tsv (Z axis @ 5s)

Harmonic series detection in Z axis

# show_harmonic_series(z_spectra, 'z', T_SEC)

show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

bvs_1_hore.tsv (Z axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.216447	1
1	73.2422	0.0392826	3
2	90.332	0.0297253	3.7
3	195.312	0.0224701	8
4	798.34	0.0184405	32.7

bvs_2_hore.tsv (Z axis @ 5s)

• Fundamental frequency: 24.4140625 Hz

	f	y	1x
0	24.4141	0.219855	1
1	58.5938	0.0495791	2.4
2	109.863	0.0368702	4.5
3	324.707	0.0300992	13.3
4	393.066	0.0230425	16.1

bvs_3_motor.tsv (Z axis @ 5s)

• Fundamental frequency: 925.29296875 Hz

	f	y	1x
0	925.293	0.0454423	1
1	292.969	0.026888	0.316623
2	1123.05	0.0182756	1.21372
3	268.555	0.00586939	0.290237
4	1225.59	0.00519174	1.32454

bvs_4_motor.tsv (Z axis @ 5s)

• Fundamental frequency: 97.65625 Hz

	f	y	1x
0	97.6562	0.044842	1
1	24.4141	0.0232632	0.25
2	1232.91	0.0146865	12.625
3	17.0898	0.00783922	0.175
4	1145.02	0.00706891	11.725

Histogram

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

Orbitals of all cross sections

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

Orbitals of 1x harmonic frequency

x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv

t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

bvs_1_hore.tsv

bvs_2_hore.tsv

bvs_3_motor.tsv

bvs_4_motor.tsv